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The effect of strongly repulsive interactions on the tunneling amplitude of hard-sphere (HS) 
bosons confined in a simple cubic (sc) optical lattice plus tight external harmonic confinement in 
continuous space is investigated. The quantum variational Monte Carlo (VMC) and the variational 
path integral Monte Carlo (VPI) techniques are used at zero temperature. The effects of the lattice 
spacing iv/k on the tunneling amplitude is also considered. The occupancies of the lattice sites as 
a function of the repulsion between the bosons are further revealed. Our chief result is, that for 
a small number of bosons (N=8) the overlap of the wave functions in neighboring wells does not 
change with an increase of the repulsive interactions and changes only minimally for a larger number 
of particles (N = 40). The tunneling amplitude rises with a reduction in the lattice spacing. In 
addition, the occupancy of the center of the trap decreases in favor of a rise in the occupancy of the 
lattice sites at the edges of the trap with increasing HS repulsion. Further, it was found that the 
energy per particle at certain optical depths is insensitive to the number of particles and variations 
in the HS diameter of the bosons. In order to support our results, we compare the VMC results 
with corresponding VPI results. 



I. INTRODUCTION 

The tunneling of bosons in optical lattices has drawn 
considerable interest in the last few years [TJ |3J |U EJ 
EJ [101 E] due to the fact the investigation of the 
tunneling amplitude in optical lattices provides insight 
into the physics of lattice bosons. This is due to its con- 
nection to superfluidity |12j and its analogy to Josephson 
tunneling in quantum devices |13j . Particularly the role 
of interparticlc interactions in determining the tunneling 
amplitude has only been given few investigations. 

In simple elementary terms, the quantum tunneling 
of a particle through a potential barrier, such as in an 
optical lattice, can occur when its total energy E is less 
than the height of the barrier Vq. The particle is de- 
scribed by a wave packet which can penetrate the barrier 
with a finite probability. As a result , an overlap between 
two wave functions at both sides of the barrier provides 
a measure for the tunneling amplitude of the particle. 
Potential barriers of this sort have been realized in quan- 
tum devices [T3] that make use of Josephson tunneling [8 
and in the recently achieved optical trapping of bosons 
|14j . When there is more than a single atom in the po- 
tential well, the mechanism of quantum tunneling is very 
much determined by the strength of the interatomic in- 
teractions. In the strongly interacting regime, correlated 
hopping occurs where atoms tunnel in pairs and 
competes with single-particle tunneling. 

Further, one of the most important properties of 
single-particle tunneling in optical lattices is that it is 
a signature of a superfluid (SF) state [15] . and if absent, 
of a Mott-insulator (MI) state [4] [7] [16]. Pair superfluid- 
ity has also been recently discussed [15]. In an MI state, 



single-particle tunneling and phase coherence are absent. 
Consequently, the particles are unable to superflow, but 
still able to hopp from one well to the other. This is 
via the correlated hopping mechanism signalled by an 
overlap of the wave functions in neighboring wells in the 
strongly interacting regime. In fact it was Foiling et al. 
[5J who showed experimentally that strong interactions 
suppress single-particle tunneling such that second-order 
correlated tunneling is then the dominant dynamical ef- 
fect. An SF state is quite the opposite, where consid- 
erable single-particle tunneling and phase coherence are 
observed. 

In this paper, we chiefly investigate the effects of 
strongly repulsive interactions on the correlated tunnel- 
ing amplitude of bosons in an inhomogeneous system con- 
fined by a tight combined harmonic optical lattice. Here, 
the external harmonic trap introduces an inhomogeneiety 
in the atomic density distributions. The tunneling am- 
plitude is measured by the overlap integral of two wave 
functions in neighboring wells obtained from the inte- 
grated optical densities. A key point in our investigation 
is to also explore the possibility of a SF to MI transi- 
tion in the presence of this inhomogeneiety since it has 
been shown that external confinement added to the op- 
tical lattice suppresses the MI state [T71[TB]. We partly 
check this case for few-boson systems. We use quantum 
variational Monte Carlo (VMC) and variational path in- 
tegral (VPI) Monte Carlo techniques in continuous space 
at zero Kelvin. To the best of our knowledge, previ- 
ous work did not consider the effects of inhomogeneiety 
on the tunneling amplitude of strongly repulsive lattice 
bosons, particularly by using Monte Carlo techniques in 
continuous space. 
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To this end, the tunneling amplitude of bosons in op- 
tical lattices has been investigated chiefly as a function 
of the optical depth (i.e., barrier height Vo) H2J [19] 
and number of particles [3 . An investigation most rele- 
vant to our work is that of Shams and Glyde [T5] . They 
evaluated the BEC density and superfluid fraction of HS 
bosons confined in an external periodic potential using 
PIMC in order to shed further light on the connection be- 
tween BEC and superfluidity. In part, they investigated 
the hopping parameter as a function of Vo, and showed 
that if Vo is increased sufficiently, the condensate is lo- 
calized into islands inside the potential wells suppressing 
superflow substantially. Further, they found that their 
external potential suppresses the superfluid fraction at 
all temperatures. 

In our investigation a key point is that, in contrast 
to Shams and Glyde, Vo is kept fixed while the repulsive 
interactions between the bosons are varied. Further, the 
effects of these interactions on atom correlations, optical 
density, occupancy of lattice sites, onsite interaction en- 
ergies, and the total energies are explored. In addition, 
the effects due to lattice spacing and number of particles 
are further revealed. 

As outlined in the method section, the bosons are 
represented by hard spheres (HS) of diameter a c whose 
repulsive interactions can be modified by changing a Cl 
thereby mimicking the Feshbach resonance technique 
[20] . Our chief result is that, for a small number of 
bosons, the overlap of wave functions in neighboring po- 
tential wells does not change with increasing HS repulsion 
and changes only minimally for N = 40. The localized 
wave functions in the potential wells do not broaden with 
an increase in a c , in contrast to the case of HS bosons in 
pure harmonic traps [5T] EH)- In the latter, the width of 
the spacial many-body wave function extends to several 
trap lengths, whereas in our case, the wavefunction in 
each well extends only slightly beyond two trap lengths 
even at large repulsion. We thus have evidence to sug- 
gest, that the optical lattice barriers prevent the wave 
functions in each well from expanding to a comparable 
extent as in pure harmonic traps. We further found that 
the energy per particle is relatively insensitive to the 
number of particles N and variations in a c as compared 
to HS bosons in pure harmonic traps [5TJ [55] , and that 
the tunneling amplitude increases with the reduction of 
the lattice spacing. In order to provide further support 
to our findings, we compare our VMC results with the 
corresponding VPI results. 

Previous theoretical work on optical lattices is abun- 
dant, particularly the SF to MI transition [7[ l^j [TB[ |P71 |531 
[Ml 013 (HI and coherent matter waves ftM |2"?1 1551 |5j?[ |3U] 
have been investigated extensively. Other investiga- 
tions included supersolidity [31], multicomponent sys- 
tems [32 , 33] , vortices [31] , solitons in a radially periodic 
lattice [3S], and p-h excitations in Mis [TB] , 

Various techniques and methods have been used to 
investigate bosons, fermions, or mixtures of them in op- 
tical lattices: the Bose-Hubbard model [3 [T7] [55] [33] as 



well as the Fermi- Hubbard model [33 [3H] in conjunction 
with Monte Carlo techniques [13 H3 H3 [33 [3J5], and 
particularly the Worm Algorithm [38] [40j [41] have been 
applied extensively. Other important techniques such as 
the Gross-Pitaevskii equation (GPE) [3J[3-S], variational 
approaches [3 1231 135] . density matrix renormalization 
group [30J , and path integral approaches [42] have also 
been used. Most of the above methods use a discrete 
space approach, whereas we evaluate the properties in 
continuous space. 



II. METHOD 

We thus consider N bosons on a combined harmonic 
optical cubic lattice (CHOCL). It consists of a 3D°^ sim- 
ple cubic (SC) lattice of Nl = 3x3x3 sites embedded 
in a tight external harmonic trap of frequency u}y, and 
trap length a; lo = y^h/muj^,, where m is the mass of the 
bosons, and T% is Planck's constant. The lattice spacing is 
given by d = ir/k, where k is the wave vector of the laser 
light. The bosons are modelled by hard spheres (HS) of 
diameter a c and their interactions Vint ( r ) ar e represented 
by a hard-core potential of diameter a c given by 

-it / \ f oo : r < a c M , 
= { : T > a c ' W 

where r is the distance between a pair of bosons. In the 
low-energy limit, the scattering between the bosons is 
purely s-wave with scattering length a s . In this limit, a s 
equals a c [5TJ [55] and the repulsive interactions between 
the bosons are modified by changing a c . Within this 
framework then, the tunneling amplitude J and the rest 
of the properties are measured as functions of a c . 

To set the stage, then, we first define the Hamilto- 
nian, and the trial wave function. Then the tunneling 
amplitude is measured by the overlap integral I over i ap of 
two wave functions in neighboring wells centered at po- 
sitions R„ and R, l+1 measured from the center of the 
trap, where n is an arbitrary integer. Here n = (pqr) are 
site indices with locations R„ = (pi + <?j + rk) in units of 
the lattice spacing d. By evaluating the overlap integral, 
the purpose is just to obtain a qualitative measure for 
the tunneling amplitude. Next the average onsite inter- 
action energy (U n ) and the average occupation number 
per lattice site (N pqr ) at position R„ = (pqr) are further 
defined. These are evaluated using VMC and VPI. We 
do not explain the VMC and VPI methods here as they 
can be found in the abundant literature [21] [43] HU [45] . 



A. Hamiltonian 

The Hamiltonian for our systems is given by 
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h 2 N N 

2—1 i—1 i<Cj 

(2) 

where Vhoi. r i) — ^ rnu! ho r i ^ s an external harmonic trap- 
ping potential with = (ccj, j/j, z\) representing a single 
particle position, m the mass of the bosons, Vi„t(ry) with 
r,j = |rj — rj| the pair-interaction potential Eq.([T]) above, 
and 



V opt (ri) = V [sm^fca:,-) + sin 2 (fc?/;) + sin 2 (fcz l )] , 

(3) 

is the optical lattice potential [7] with V$ the optical 
depth. Essentially, V op t(vi = R„) = at the lattice-site 
positions R„ which are the locations of the potential-well 
minima of Eq.j3|. That is, V op t(Rn) = Vb[sin 2 (p7r) + 
sin 2 ((77r) + sin 2 (r7r)] = where k • R„ = (integer) x n. 
Thus the lattice-site positions are implied in V op t(ri) 
by its very construction. Experimentally, the optical 
lattice potential is obtained from a superposition of 
three pairs of mutually perpendicular, counterpropagat- 
ing laser beams of intensity proportional to Vq- For a SC 
lattice, Vq must be the same for each direction. In this 
paper, we write energy and length in units of the trap, 
huJho and ciho = \Z^/fn^ho, respectively. 



B. Trial wave function 

The many-body trial wave function is given by 

N 

*({r},{R}) = nex P (-«r 2 )^(r 2 ,{R}) [] / (l^l), 

(4) 

with t/j(ri, {R}) a Wannicr-likc function defined as 

N L 

<A(r,,{R}) = ^0( ri ,R n ), (5) 

n=0 

where {r} = {ri,r 2 , • • • , r^} is a set of spatial vectors 
describing the positions of the bosons from the center of 
the trap, {R} = {Ri, R 2 , • • • , Rjv l } is a set of vectors 
describing the positions of the lattice potential minima 
on the 3x3x3 cubic lattice cage considered, and a is 
a variational parameter signalling the strength of the ex- 
ternal harmonic confinement. Essentially, a is the inverse 
overall width of the total wave function of the system in 
the CHOCL. Further, a controls the volume of the exter- 
nal harmonic confinement and, therefore, the number of 
lattice sites to be occupied away from the center of the 
trap. Eq.(|5| is constructed in a manner similar to that 
of Jin et al. |33j in that we sum over a number of local- 
ized single-particle wave functions 0(rj,R n ) centered in 



the optical lattice wells, each at R„. 0(ri,R„) is con- 
structed similarly to the wave function used by Li et al. 
[7] and is given by: 

^(r^Rn) = exp[-/3(r 4 - R n ) 2 ] x 

[1 + - X n f - a{x t - X n ) 4 ] x 

[1 + 7 (W - Yn? - a{ Vl - Y n ) 4 ] x 

[l + ^z,- Z n ) 2 -a{ Zl - Z n f] , (6) 

where /3, 7, and a are further variational parameters in 
addition to a. The onsite repulsion is partly controlled by 
the local density |<^)(ri,R„)| 2 via the variational parame- 
ter p in Eq. (|6| , which confines the particles at each lat- 
tice site (n). The interactions in the trial wave function 
are taken into account by the usual HS Jastrow function 

|ZHE2|: 

*-•«>-{ 1- £° ; w 

where r t j is the distance between a pair of bosons. 

The VMC wave function Q is optimized by mini- 
mizing the average energy with Powell's technique [46] 
employed in a previous publication [47 . The ground 
state configuration is achieved when all the particles are 
distributed symmetrically in the spherical coordination 
shells around the center of the trap, and the overall 
boson-boson repulsion is minimal. Once the trial func- 
tion Eq. Q has been optimized, it is plugged into the VPI 
code as a starting wave function. The same properties 
are then evaluated as those using VMC. 

C. Tunneling 

Since we consider in this work only the strongly in- 
teracting regime, single-particle tunneling is very much 
suppressed in our systems. Hence, the atom-pair tun- 
neling is dominant and is included in the overlap of the 
wave functions in neighboring wells. A qualitative mea- 
sure for the overlap of the wave functions is evaluated, 
for simplicity, from trial-density functions fitted to slices 
of the integrated optical densities along one axis. For 
this purpose, the trial-density functions are similar to 
the one-dimensional version of Eq. (|6j) with an additional 
amplitude factor A n , that is 

(j) flt (x,X n ) 2 = A n exp[-/3„0 - X n ) 2 } x 
\l+ lrL {x-X n ) 2 + a n {x-X n f\. (8) 

The absolute value is considered to make sure we do not 
run into negative densities. An example of such a fit is 
shown in Fig. [T] below. Here, the fitting function is: 

3 

F(x;X 1 ,X 2 ,X 3 ) = J2<( ) fit{x,Xn) 2 (9) 

n=l 
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FIG. 1: Fits using the forms Q and (jgj to ID VMC and VPI 
integrated optical density slices along the x-axis for a system 
of N — 8, Vb = 10, and k = 1.2-7T. The sum of squares for 
both fits is 3.41 x 10~ 4 for VMC and 7.41 x 10~ 4 for VPI, 
respectively. 



having peaks centered at three potential minima X\ 
(left), Xi (center), and X$ (right) and with three sets 
of fitting parameters A n , f3 n , j n , a n , and the latter X n , 
where n runs from 1 to 3. In fact, the central position 
Xi always turns out to be exactly zero as required. The 
sum of squares 



(\4>f it (xi,x n )\ 2 - \4> MC {xi,x n )Y 

•(io) 

is minimized with respect to the above sets of fitting 
parameters. Here \ 2 is a sum over all P data points 
of the one-dimensional integrated MC optical density 
iJ2n=i\ ( t ) Mc(xi,X n )\ 2 }. By this minimization, we get 
values of x 2 < 10~ 3 indicating a good fit. After opti- 
mization, the overlap integral 



overlap 



(f>fit(x, Xi) <f>fi t (x, X 2 )dx + 
<Pfit(x,X 2 )(j>fit(x,X 3 ,)dx, (11) 



is then evaluated using a simple elementary numerical 
technique. 

We must emphasize, that it is hard to describe the 
single-particle tunneling amplitude for strongly interact- 
ing systems using the well-known exchange integral [7J , 



J n = / d 3 rc/> n (r) 



2m 



V(r) 



^ + l(r), (12) 



where <j) n and <p n +\ are single-particle wave functions at 
sites n and n + 1, and V(r) a single-particle potential. 



This is because the interactions as described by the Jas- 
trow factor [Eq.d7|] are not included in Eq.(12l. In fact, 



it is anticipated that the single-particle wave function 
narrows as the HS repulsion rises, counteracting the ef- 
fects of the broadening due to the Jastrow factor as the 
number of particles is increased beyond a certain limit. 



D. Bosons-optical lattice overlap 

A measure for the extent of the overlap between 
</>(ri,R n ) and V op t(ri) for all particles i = 1 to iV and 
all lattice sites n = 1 to Nl is given by 



l boson-ol) — 



E 

71=1 



(MC) 



' 0(r,R„)y opt (r)d 3 r, 
(13) 

where J^ MC \ stands for a MC configurational integral de- 



fined in Sec II H below, and OL stands for optical lattice. 
Further (Iboson-ol) is an indirect measure of the total 
tunneling amplitude of the system into the optical lattice 
potential barriers. 



E. Discrete onsite interaction energies 

In order to explore the single-particle response to the 
change in the HS repulsions, we evaluate the onsite in- 
teraction energies at individual lattice sites. The onsite 
interaction energy for a single particle at each lattice site 
(n) is evaluated using [Jj |48] 



(Un) 



9 



-4ar 



\ 4 d 3 7 



(14) 



(MC) 



where g = ^irTi 2 a c /m with a s replaced by a c . Here the 
effects of the external trap are included via the factor 
[exp(— ar 2 )] 4 . 



F. Occupancy of Lattice Sites 

Another single-particle response is the atom-number 
occupancy of each lattice site. The average number of 
particles occupying a particular lattice site (n) is evalu- 
ated via 



N 



-2ar 



|0(r,Rr, 



| 2 d 3 r 



(15) 



(MC) 



An alternative way to determine the occupancy of each 
lattice site, is to divide the confining volume into equiva- 
lent cubic bins whose corners are sets of four lattice sites 
R n and count the number of particles collected in each 
bin. We used this counting method to evaluate the VPI 
(N n ), whereas (fl5j» to obtain the VMC (N n ). 
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G. Correlations between the bosons 

The pair correlation function g(r) is evaluated using 
VPI by binning pairs of particles in each r. The idea 
is to provide a measure for the strength of the boson- 
boson correlations and its dependency on the repulsive 
interactions in a tightly confining environment . 

H. Numerics 



The integrals (131, (14 1 and (151 are evaluated using 



standard Monte Carlo integration |45 



1 N 1 

-Y- , 

N ^ M ^ 

i=l c=l 



M 



(MC) 



0{v cl ) 



-lur' 1 . 



where c is a configurational index, i a particle index, 
and M is the total number of VMC configurations. We 
further sum over all N particles and divide by N to get 
the average. The denominator is a weight we used for 
all MC integrals. In VPI we apply the same integration 
technique, except with M taken as the total number of 
time slices used in the path integral. Each time slice 
contains one configuration of TV particle positions. 



I. Computational complexities 

As much as VPI is an accurate method for the present 
determination of the properties of lattice bosons, using 
a number of particles exceeding N — 8, this method be- 
gins to constitute a heavy-computational and CPU-time 
consuming technique. Particularly, for large optical po- 
tential barriers Vq > 10 fiujho, the bosons take a very 
long CPU time to tunnel through the barriers in order 
to achieve a symmetric optical density distribution inside 
the whole trap. We therefore limited the number of par- 
ticles to N = 8 when using VPI. On the other hand, the 
less accurate VMC method allows the bosons to diffuse 
quickly through the barriers providing a computationally 
cheap and fast technique to investigate qualitatively the 
properties of bosons in optical lattices. Our aim is, how- 
ever, not to present methods but to concentrate on the 
physics of quantum tunneling. 



III. RESULTS 

In what follows, we present the results of our calcula- 
tions. We particularly focus on the effect of interactions 



on the tunneling amplitude cx /, 



overlap- 



We further ex- 



plore the effect of a c on the onsite interaction energies 
(U n ), energy per particle (E)/N, the occupancy of lattice 
sites (N n ), the correlations g(r), and the optical density. 
As mentioned before, we compare our VMC results with 
corresponding VPI results. 
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,FIG. 2: Tunneling amplitude oc I over i av [Eq.( 11 1] vs the hard- 
sphere (HS) diameter a c , obtained from fits to VPI optical 
densities. The system is a HS Bose gas of JV = 8 particles in 
a 3 x 3 x 3 cubic optical lattice of depth Vb = 10 embedded 
in a tight external harmonic trap of trapping frequency loho- 
Open circles: k = tt, open triangles: k = 1.27T, and open 
diamonds: k = 1.4-7T. Energies and lengths are in units of the 
trap. 
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FIG. 3: As in Fig. [2] but using VMC for TV = 40, V = 20 
and the indicated values of k. 



A. Tunneling 

loveriap is now displayed as a function of a c . Fig. 
[2] displays the VPI results for systems of N = 8 and 
Vq = 10 with k = 7T (open circles), k = 1.2-k (open tri- 
angles) and k — lAir (open diamonds). For these values 
of N and Vq , I overlap shows a very weak response to the 
effects of increasing a c , since it seems to remain practi- 
cally constant. As k is increased, loveriap rises indicating 
a profound effect of the lattice spacing d — ir/k on the 
tunneling amplitude. 

The case for N = 40 and Vb = 20 is shown in Fig. 
[3] for A;— values of tt and 1.27T (solid circles and triangles, 
respectively). For k — 1.2w, loveriap decreases overall 
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FIG. 4: VMC parameter (3 of the trial wave function [Eq.pj)] 
as a function of the hard-sphere diameter a c for two HS Bose 
gas systems in a combined harmonic simple cubic optical lat- 
tice of 3 x 3 x 3 sites. Solid circles: HS Bose gas of Fig. [2] at 
k — 1.27V. Open triangles: HS Bose gas of Fig. [3] at k = tt. 



FIG. 6: As in Fig. [5] but for 
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B. Width of the single-particle wave function in 
each potential well 
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FIG. 5: Average VMC onsite repulsive interaction energy 
(U n ) at each lattice site n for a HS Bose gas of N = 8 parti- 
cles, optical lattice depth Vb, and HS diameter a c — 0.02 in 
a 3 x 3 x 3 cubic optical lattice embedded in a tight exter- 
nal harmonic trap of trapping frequency uiho- Energies and 
lengths are in units of the trap. 



with increasing a c , whereas for k = I. On it rises some- 
what up to a c — 0.18 and then begins to drop. However, 
the response is still weak and insignificant. Hence, a 
slightly stronger response is obtained for the tunneling 
amplitude when using a larger N. Further, the effect of 
changing d = ir/k on I overlap is much more pronounced 
than changing a c . Signals for a stronger response of our 
systems to the rise of a c are found in the single-particle 
properties such as the width of the single-particle wave 
function and the onsite interaction energy. 



As already indicated in the Method section, the width 
of the single-particle wave function </>(r,R„) [Eq.Q] in 
each well, is expected to decrease with a rise of the HS re- 
pulsion in order to counteract the effects of a broadening 
due to the Jastrow function for a larger number of parti- 
cles. The width, or better the inverse of it, is described 
by the VMC parameter (3 in Eq.([6]). In Fig. [3] the VMC 
parameter [3 of the local wave function, 0(r,R„), is dis- 
played as a function of a c . The open triangles display (3 
for the HS Bose gas of Fig. [2] at k — \.2ti\ solid circles: 
HS Bose gas of Fig. [3] at k = tt. It is found that this 
width (cx l/\/]3) decreases with the rise in the repulsion 
for N = 40 and remains almost constant for N = 8. 



C. Onsite interaction energies 

Discrete distributions of (U n ) over all 27 lattice sites, 
for systems of N — 8, Vb = 10, k = tt, and two cases 
a c = 0.02 and 0.10, are displayed in Figs. [5] and [6] re- 
spectively. The abcissa display the site index n. All 
(U n ) which have the same magnitude within a small mar- 
gin of error, belong to the same set of nearest neighbors 
from the central lattice site. That is, moving vertically 
from the tip of the discrete distribution down to the bot- 
tom, we move from the center to the first, second, and 
third set of nearest neighbors, respectively. The distri- 
bution is uniform around the trap center. The effects 
of the external trap are manifested in the magnitude of 
(U n ) which is highest at the center and declines towards 
the edges of the trap due to the factor exp(— 4cvr 2 ) in 
Eq.(|i"4"l). Obviously, a large percentage of the particles 
is concentrated at the central lattice site. (It is antici- 
pated then, that if the strength of the external harmonic 
trap is increased further, this causes a drop in the density 
~ exp(— 2or 2 ) |</>(r,R n )| 2 at the lattice sites towards the 
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FIG. 7: Average onsite repulsive interaction energy {U( pqr )) 
vs the hard-sphere (HS) diameter a c at four lattice sites rep- 
resentative of the whole lattice. The system is the HS Bose 
gas mentioned in Fig. [2] at k — n. Solid and open circles: 
VMC and VPI results for (f/(ooo)}- Solid and open triangles: 
likewise for (00 — 1). Solid and open diamonds: for (Oil). 
Solid and open inverted triangles: for (111). 



edges of the trap and a certain rise at the center of the 
trap.) Eventually, all the bosons will occupy the central 
lattice site and only if their HS diameter is of a magnitude 
that allows all of them them to fit inside one lattice-site 
volume ~ 4-7r(<i/2) 3 /3. Thus, a large amount of repul- 
sive potential energy would be stored in the trap center. 
The patterns of Figs. [5] and [6] show also that the bosons 
minimize their potential energy arising from the external 
harmonic trap by maximizing their occupancy of the trap 
center, and minimizing it towards the edges of the trap. 

In addition, (U n ) is displayed as a function of a c in 
Figs. [7] and [8] for the systems of Fig. [2] at k — ir and 
1.27T, respectively, and for four lattice sites: R„ = (000), 
(00 — 1), (111), and (Oil), representative of the whole 
lattice. In general, {Utp qr \) rises with a c signaling a rise 
in the repulsive energy of the particles in each well. The 
rate of growth of (C/( P9r )) is largest for the trap center 
and declines towards the edges of the trap. 



D. Boson-optical lattice overlap 



FIG. 8: As in Fig.[7] but for k = 1.2tt. 



In this section, we evaluate (Ibosons-ol) [Eq. (13 1] 



which measures the overlap of all the single particle wave 
functions in the wells with the periodic optical lattice po- 
tential. In essence, it measures also the total tunneling 
amplitude of the particles into the optical lattice poten- 
tial barriers. Figs . [9] and [To| display (Iboson-ol) vs. a c 
for the systems of Figs. [2] and [3] respectively. One can 
see that (Iboson-ol) is practically invariant for N = 8, 
but the system with TV = 40 reveals a slightly stronger 
response to changes in a c . For example for N = 8 and 
k = 7T, the change in (Iboson-ol) from a c = 0.02 to 
0.3 is only by —1.45%. The other changes for N — 8 are 
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FIG. 9: VMC disorder integral [Eq.(13l] vs the hard-sphere 
diameter a c for the systems in Fig. [2 
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FIG. 10: VMC disorder integral [Eq.Q] vs the hard-sphere 
diameter a c for the systems in Fig. [3] 



of the same order of magnitude. This in turn indicates, 
that the width of each localized wave function for these 



systems is practically invariant with the rise in repulsion. 
This may come as a surprise, since it is known [2TJ [22] 
that at very large HS diameters such as the ones used 
here, the wave function ought to become very broad in- 
deed. It seems that the localization effect of the optical 
lattice potential overwhelms the repulsive forces tending 
to broaden the wave function in each well. 

One further observes that (Iboson-ol) increases 
with k in agreement with the results of Figs. [2] and [3] 
where I overlap rises with increasing k. 

Further, we note that in Fig. 10 (Iboson-ol) rises 
notably with a c indicating a rise in the total tunneling 
amplitude into the optical potential barriers. We were 
unable to obtain values of (Iboson-ol) for k = 1.2tt be- 
yond a c = 0.18 as we could not find VMC energy minima 
for these systems. We must emphasize that the overlap 
of the bosons with the optical lattice potential is a dif- 
ferent measure than I over i ap , which measures the overlap 
of two wave functions in neighboring wells. 



E. Occupancy of lattice sites 

Figure [TT] displays the average occupancy of four lat- 
tice sites (N(pqr)) for systems of Fig. [2| (pqr) = (000), 
(00 - 1), (011), and (111). Each frame is labelled by the 
corresponding lattice site. The solid and open diamonds 
display VMC and VPI results for k — 1.4tt, whereas the 
solid and open triangles those for k = 1.2ir, respectively. 
The VMC results have been obtained by Eq. ( fl5| , whereas 
the VPI results from simple counting of the bosons in cu- 
bic bins. In frames (000) and (00 — 1), the occupancy 
decreases with a c , whereas in frames (011) and (111) it 
rises. Further, the occupancy is higher for a lower fc, i.e., 
a larger lattice volume. The VMC and VPI results almost 
match for (000) and (00 — 1) and are close to each other 
for (011) and (111). There is good agreement between 
the occupancies calculated by Eq.( 15 1 and those obtained 



by counting. This indicates that the single-particle wave 
function in Eq.([6| is suitable to describe the bosons at 
each lattice site. 



F. Correlations 

The VPI pair correlation functions g(r) for the sys- 
tems of N — 8, Vo = 10, k — it, and various a c in the 
range 0.02 < a c < 0.30 are displayed in Fig. 12 It is 



observed that the correlations are strongest for a c = 0.02 
and gradually weaken as a c is increased. 
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FIG. 11: (Color online) Average occupancy {Nr pqr )} vs the 
hard-sphere diameter a c at the four lattice sites (pqr) indi- 
cated, which are representative of the whole optical lattice. 
The system is the HS Bose gas of Fig. [2] Solid and open 
diamonds: VMC and VPI results for k — IAtt, respectively. 
Solid and open triangles: likewise, but for k = 1.2n. From top 
to bottom frame: R13 = (000) is the center, R12 = (00 — 1) 
the first, R17 = (011) the second, and R26 = (HI) the third 
nearest neighbor to the center, respectively. Site indices n are 
according to Fig. [5] 



G. Energy 

The average Monte Carlo (MC) energies per parti- 
cle (E)/N as functions of a c for the systems of Figs. [2] 
and [3] are displayed in Figs.[l3j and [14] respectively. 
Additionally, the energies for a system of N = 2 and 



Vo = 2, are displayed in Fig. [15] as well. The same 
legends are used in the latter three figures as follows. 
Solid and open diamonds: VMC and VPI results for 
k = IAtt, respectively. Similarly for triangles and cir- 
cles: k = 1.2-k and n, respectively. We note that, for 
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FIG. 12: VPI pair correlation function g(r) for the indicated 
a c . The system is the HS Bose gas of Fig. [2] at k — tt. 
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FIG. 14: As in Fig. [13] but for the HS Bose gases of Fig. [3] 
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FIG. 13: Average energy per particle (E) /N vs. a c for a 
HS Bose gas system of N = 8, Vo = 10, and the indicated 
values of fc. Solid and open diamonds: VMC and VPI results 
for k — IAtv. Solid and open triangles: likewise, but for 
k — 1.27V, and solid and open circles: k = tt. 



most of the systems considered, (E)/N changes almost 
linearly within the given range of a c at a much slower 
rate than HS Bose gases in pure harmonic traps [21] [22] . 
For example for TV = 2 with k = tt, the VPI (E)/N 
rises only by AE VPI /N - +2.98% from a c = 0.02 to 
a c = 0.3, despite the fact that there is a very large change 
in a c (by a factor of 14!). Similarly for N — 8 and 
k = tt, AE VPI /N is - +4.48%. However, AE VPI /N 
for k — \Att is more pronounced, being ~ +6.26%. The 
corresponding VMC results are as follows: for N = 2 

with k = tt, AE VMC /N h4.48%, for N = 8 with 

k = n, AE V Mc/N h5.20%, and for N = 8 with 

k = 1.4tt, AE VM c/N 1-12.37%, respectively. For 

N = 40, AE/N changes at a much faster rate than for 
a lower N. Further, the energy (Evmc)/N rises with 
increasing k. 
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FIG. 15: As in Fig. [TJ but for iV = 2, Vo = 2, and the 
indicated values of k. 



H. Optical density 

In this section, the integrated VPI optical density 
(OD) is displayed at different a c . The goal is to explore 
the effects of variations in the interactions on the optical 
density profiles and on the tunneling amplitude. Par- 
ticularly we show visually that, for a small number of 
particles, the width of the wave function in each well 
does not broaden with an increase in a c even up to order 
0.1 and that the overlap between the wave functions does 
not change either. 

Fig. [Tg] displays the integrated VPI OD for N = 2, 
Vo = 2, and k — tt. From the top to bottom frames, 
a c = 0.02, 0.14, and 0.30, respectively. As expected, the 
amplitude of the density is maximal at the center of the 
trap and lowest for the lattice sites near the edges of 
the trap. Further, the tunneling of the system into the 
potential barrier of the external harmonic trap is sup- 
pressed due to the steep rise of this barrier as one goes 
away from the center of the trap. Since the optical depth 
is very shallow, the wave function of the system is spread 
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FIG. 16: Integrated VPI optical density for the HS Bose gas 
system of N = 2, Vb = 2, and k = n, at three values of 
a c . From top to bottom frames: a c — 0.02, 0.14, and 0.30, 
respectively. Density and lengths are in trap units. 



over the whole lattice, indicating the dominance of the 
tunneling process [1] over the localization effects of the 
optical lattice potential wells. Superfiow in these systems 
is thus prevailent. 



Fig. 17 displays one-dimensional density slices of Fig. 
16 along the x-axis, and additionally three functions of 
the form ^ fitted to this density profile. The fits reveal 
the extent of the overlap between the wave functions in 
neighboring wells. It is observed, that this overlap does 
not change with a rise in a c . The left (dashed line) and 
right (dotted line) density functions are much broader 
than the central function (thin dashed line) and all three 
peaks have almost the same amplitude. In fact, the over- 
lap here is between three wave functions as those at the 
edges are able to spread over several trap lengths. 

Fig. [18] is the same as Fig. [16] but for N — 8 and 



FIG. 17: Integrated one-dimensional optical density slices 
of Fig.[l6] along the x-axis. Solid line: fitting function 
f(x) = F(x; Xi,X2,X 3 ) Eq.{9|, open circles: VPI data, thick 
dashed line: left fit \4>fit(x,Xi)\ 2 , thin dashed line: central fit 
\4>fit{x,X 2 = 0)j 2 , dotted line: right fit \4>j it (x, X 2 ) \ 2 ■ From 
top to bottom, same a c at in Fig. |16| 



Vb = 10- Here, the bosons are more localized at the 
lattice sites, and the tunneling process is less dominant 
than the localization effects. Again, the effects of the 
external trap are manifested: the density diminishes to- 
wards the edges of the trap. Fig. [19] displays again the 
one-dimensional slices of Fig. [T8] along the x-axis. Here, 
the overlap is much less pronounced than in Fig. |17| and 
the amplitude of the central peak is larger than the peaks 
at the edges. The left and right peaks do not overlap as 
in Fig. [17] It seems that a slight increase in the num- 
ber of particles has a profound effect on the density pro- 
files in that it changes the ratio of the amplitude of the 
central density to that of the density at the trap edges. 
Further, the width of the densities at the edges of the 
trap for N — 8 has dropped by almost 2 trap lengths as 
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compared to N = 2. Nevertheless, there is still overlap 
between the wave functions in the neighboring wells even 
at high repulsion. For N = 8, the amplitude of the cen- 
tral peak declines with increasing a c . The change in the 
height of the central peak as one increases a c from 0.02 to 
0.14 is ~ 33.3%. Remarkably, even with strongly repul- 
sive bosons, the density of the system in each well largely 
peaks at the potential minimum of each well instead of 
being broader and more evenly distributed around each 
lattice site as in Fig. [16] 

The latter situation does not change very much for a 
larger number of particles. Fig. 20 displays the integrated 
optical density for the system of N — 40, Vq — 20, and 
k = l.2w and Fig. [2T]its density slices, as in the latter fig- 
ures for the optical density. Again, the peaks in each well 
do not broaden very much with a rise in a c . The central 
peak looses amplitude in favor of the peaks at the edges. 
On closely inspecting the area of the overlap between 
neighboring peaks, one can depict a small decline in the 



tunneling amplitude of the system as a c is increased. In 
Figs. [19] and |2T] the positions of the left and right peak 
maxima shift to the left and right, respectively, as a c is 
increased. 



I. Alternative VMC Trial Function 

In order to test the validity of using a Gaussian 
(exp(— ar 2 )) in our trial wave function Eq.Q, we cal- 
culated the optical densities for N — 8, Vq = 10, and 
k = 7r using a slightly modified version of our trial wave 
function: 



*({r},{R}) 



N 



n exp (~ 



er?)V(ri,{R})n/(ki- r jl). 
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with an additional parameter e added to the previous set 
of variational parameters. Nevertheless, this produces 
similar results to the previous section. Fig. 22 displays 
the density slices for the system of N = 8, Vb = 10, and 
k = 7r, using Eq.(17l. The wave functions in each well 



FIG. 21: As in Fig.|r7]but for TV = 40, V = 20, and k = 1.2tt 



do not broaden, but the central peak looses amplitude 
with increasing a c , whereas the peaks at the edges gain 
amplitude. 



IV. SF TO MI PSEUDOTRANSITION 

Although our chief goal was to fix Vq and vary a c only, 
we nevertheless divert a little and explore the evolution 
of a HS Bose gas in a CHOCL by changing Vb and fixing 
a c . The goal is to check whether there is a critical Vq at 
which a SF to MI transition occurs in the presence of an 
external harmonic trap. 

Fig. [23] displays the evolution of the integrated opti- 



cal density (map view) for a system of N — 8 particles, 
k = 7r, fixed HS diameter a c — 0.14, and various optical 
depths Vb- From top to bottom: Vq = 1, 2, 6, and 10, 
respectively. As Vq is increased, the overlap of the wave 
functions in neighboring wells drops, however, the change 
is not abrupt but gradual. That is, there does not exist 
a critical point for a superfluid to Mott-insulator transi- 
tion. 



again 



Fig. 24 displays, similarly to Figs. 17 and 19 
the density slices of Fig. [23] along the x-axis, where one 
can clearly see how the overlap, and therefore the tun- 
neling, drops. We anticipate then, that as Vb is increased 
further, the tunneling will drop substantially. Note, how- 
ever, that the amplitude of the central peak rises with 
increasing Vq. 
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FIG. 22: As in Fig. |19| but using the alternative Gaussian in 
Eq.fTFf (Alt. Gauss). 



DISCUSSION AND CONCLUSION 



In summary, we have investigated the effects of 
strongly repulsive interactions on the tunneling ampli- 
tude of HS bosons confined in a 3D a ^ simple cubic optical 
lattice plus by a tight spherical harmonic trap. The tun- 
neling amplitude was measured by the overlap integral 
of the wave functions in neighboring wells. The external 
harmonic trap was introduced in order to cause an in- 
homogeneiety in the density distribution of the trapped 
system. Another goal was also to explore the effect of this 
inhomogeneiety on the SF to MI transition of bosons in 
optical lattices. 

It was found that for small N, the tunneling amplitude 
does not change by increasing the boson-boson repulsion 
at fixed lattice spacing d = Tt/k. This is mainly be- 
cause the width of each localized wave function does not 
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FIG. 23: Map view of integrated VPI optical density for vari- 
ous optical depths Vq. The system is a HS Bose gas of iV = 8, 
fixed a c — 0.14, and k — l.On which is confined in a CHOCL. 
From top to bottom: Vo = 1, 2, 6, and 10. Density and 
lengths are in trap units. 



broaden with a rise in the repulsive forces. The tunneling 
amplitude, however, rises by reduction of d at fixed a c . 
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FIG. 24: Density slices of Fig. [23] along the x-axis. The leg- 
ends describe the same functions as in Fig. |17| 



Even at very large repulsion, tunneling is still present in 
the system and does not vanish completely. The pres- 
ence of an external trap seems to suppress the SF to MI 
abrupt transition. The transition is gradual and there is 
no criticality. 

The energies change only weakly with increasing a c 
for N — 2, but a more pronounced change is observed for 
N = 8, and is much larger for N — 40. The occupancy 
of the central and first nearest neighbor sites drops with 
increasing a c , and rises at the second and third nearest 



neighbor sites. 



A. Occupancy and onsite interaction energies 

Although the occupancy of particles in Fig. [TTJat sites 
(000) and (00 — 1) drops with increasing a c causing a 
drop in the density at the central lattice site, the corre- 
sponding (C/(ooo)) an d (J7(oo-i)) still r i se - Thus the effect 
of increasing g via a c on (U n ) overwhelms that of the 
drop in the occupancy at each lattice site. The net re- 
sult is that (U( pqr )} increases with increasing g. Further, 
the occupancy of the lattice sites plays a crucial role in 
determining the relative magnitudes of (C/( p9r )). Since 
(000) has the highest occupancy at all a c , in Figs. [7] and 
[8] it has also the highest onsite repulsion. Then with de- 
creasing occupancy, the onsite repulsion decreases in the 
order t/ (00 o), f7( 00 -i), U (ml) , and U (111) at all a c . In Fig. 
[IT) one can also see that (000) has the highest, whereas 
(Til) the lowest occupancy. 

In equilibrium, one could assume that there are as 
many particles tunneling into a lattice site as out which 
preserves the occupancy at each lattice site in equilib- 
rium 49J. When a c is changed, the particles relocate 
themselves in order to reach a new equilibrium configu- 
ration. If a c is increased, particles are pushed out of the 
central site and occupy others until a new equilibrium 
configuration is reached again. Due to the external trap 
they are not able to tunnel further away from the edges 
of the trap. Further, note it is possible that during the 
equilibration process the particles having tunneled from 
the center of the trap to the first neighbors continue tun- 
neling to the second and third. In fact, it would be in- 
teresting to know how the particles' motion is channelled 
along the sites of a larger simple cubic optical lattice, as 
the repulsion between the bosons is increased. Finally, 
the average (noninteger) occupation numbers (N^ pqr ^) at 
each lattice signal that the systems are SF. This indicates 
also that the particles are fragmented all throughout the 
lattice. 



B. Optical density and tunneling 

As we have seen in Figs. [16]- [21] as much as the re- 
pulsion was increased, the overlap between neighboring 
wells did not change significantly, even at large N. But 
the change in the overlap is more pronounced if one in- 
creases Vq while keeping a c fixed. This was the situa- 
tion displayed in Fig. [24] We conclude then, that a MI 
state can be obtained more efficiently by rather increas- 
ing the optical depth Vq than the HS repulsion between 
the bosons, a c . 

Similarly to our findings in Sec IIIH Greiner et al. [1 



showed, by applying the Bose-Hubbard (BH) model to a 
3D optical lattice, that if the tunneling dominates the 
BH Hamiltonian, the single particle wave functions are 
spread over the whole lattice, and phase coherence ex- 
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isting between the lattice sites forms a SF. On the other 
hand, if the atom-atom interactions dominate, there is 
no phase coherence. Therefore, the single-particle wave- 
functions become localized at their lattice sites with a 
fixed number of atoms, thus forming a MI state. In our 
case, however, we do not get a pure MI even at very high 
repulsion because tunneling is still present there. In the 
highly repulsive regime, we may nevertheless be talking 
about a quasi MI state where tunneling is present but 
vanishingly small. 

In Fig. |23[ if one keeps increasing Vq, we anticipate 
that the tunneling amplitude will eventually drop to zero 
signaling initially what looks like an entrance into the MI 
regime. However, in order to identify a pure MI regime, 
integer occupancy of each lattice site is needed. Because 
the number of particles is less than the number of lattices 
sites (N/Nl < 1), the particles are fragmented [25] and 
we do not get integer occupancy at each lattice site. This 
was shown in Scc |IIIE] and as a result this system may 
be regarded a quasi MI with a small SF component. 



due to the presence of external confinement. This forces 
an inhomogcncicty in particle distribution over the lattice 
sites. 



E. Energy 



The rise in (E/N) with k in Figs 13 



14 



and 



15 



could be attributed to Heisenberg's uncertainty princi- 
ple. A decline in the available volume the bosons can oc- 
cupy in the optical lattice wells reduces the uncertainty in 
their positions. Thus, the uncertainty in their momenta 
increases and as a result their kinetic energy (quantum 
pressure) rises. It can also be seen that a change in k has 
a more profound effect on E than a change in a c . There- 
fore, the lattice dimension seems to play a profound role 
in determining the SF properties of bosons in optical lat- 
tices in that it can substantially control the tunneling 
amplitude between neighboring wells. 



C. The role of the lattice spacing d — ir/k 

In Figs. [2] and |3] the reduction in n/k reduces the 
width of the optical potential barriers, thus reducing the 
distance a particle needs in order to tunnel completely 
through a barrier. As a result, an increase in the overlap 
occurs between two neighboring single-particle wavefunc- 
tions </>(r,R n ) and <^(r,R n +i) causing a rise in Ioveriap- 
The increase in I over iap with k could be further explained 
on the basis of a recent evaluation of the tunneling rate 
in ID through a single, tilted optical barrier without ex- 
ternal confinement by Huhtamaki et al. [3 . If a and b 
are the left-hand and right-hand classical turning points 
of a potential barrier, then if the distance between a and 
b is reduced by increasing the lattice wave vector k, the 
overlap between the wave functions in neighboring wells 
rises. 



D. The effects of the external harmonic trap 

The mobility of the bosons is eventually restricted to 
the confining volume of the external harmonic trap. By 
increasing the size of the bosons with a c , the available 
free space for motion in the confined volume is reduced. 
Further, as the repulsion rises with a c , the tunneling am- 
plitude remains constant for a small number of particles 
as they get locked inside the lattice sites. A tight 
external harmonic potential also suppresses the critical 

In 



SF to pure MI transition as was shown in Fig. 24 



fact, Gygi et al. [T7] argued similarly that the presence 
of a gradient due to the external confining potential de- 
stroys quantum criticality in the transition from the SF 
to the MI phase; this gradient is present in our simu- 
lations. Further, as high as the repulsion between the 
bosons in our systems becomes, tunneling is still present 



F. Other work 

Similarly to Bach and Rzazewski |24j . we used Gaus- 
sians of the form exp[— a(x — x n ) 2 ] but weighted by a 
polynomial as in Eq. ^ . These authors determined that 
the mobility of the atoms, as described by the hopping 
parameter t and their interactions U , is controlled by the 
optical-lattice depth. If this depth is shallow, the atoms 
become delocalized over the whole lattice as it almost 
occurs in Fig. [16] 

Li et al. [7 investigated the SF to MI transitions in 
atomic BECs confined in optical lattices by using the 
BH model. Using an isotropic cubic lattice, they have 
chosen a variational trial function of the form g[u) — 
(1 + cm 2 ) exp(— (3u 2 ) (a Wannier function); our single- 
particle wave function is the same as theirs, except for an 
additional factor cx u 4 [Eq.(|6])]. They addressed the pos- 
sibility of observing SF to MI transitions for an average 
lattice-site occupancy larger than one. It was noted that 
by increasing the repulsion between the lattice bosons, 
their wavefunction in each well broadens, thus enhanc- 
ing J between neighbouring lattice sites. In contrast, 
our atomic clouds did not expand due to the presence of 
external confinement. Li et al. evaluated the on-site in- 
teraction energy by variationally minimizing the energy 
with respect to the parameters a and f3. 

Capello et al. [23] used VMC to describe the MI tran- 
sition (MIT) by means of a variational Gutzwiller wave- 
function in discrete space and optimized the variational 
parameters to achieve the ground state system. In con- 
trast we used a wave function in continuous space and 
the total energy was minimized in our research as they 
did. 
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